## albma: A Package for Analysis with Atheoretic Lags
## Skyler Cranmer and Doug Rice
## This package implements the methods described in
## Skyler J. Cranmer, Doug Rice, and Randolph M. Siverson "What to do about Atheoretic Lags"



## For getting CIs 
output.summary <- function(mod){
	mod.table <- cbind(summary(mod)$statistics[,1:2],
	(summary(mod)$quantiles[,1]),
	(summary(mod)$quantiles[,5]))
	dimnames(mod.table)[[2]] <- c("Mean","Std. Dev.","2.5%","97.5%")
	mod.table }


## Lead Function
lead <- function(var, lead, by, data, newvar=FALSE, newname){
	if (missing(data)) 
		data <- environment(formula)
	if(newvar == TRUE && missing(data))
		stop("A data frame must be specified if newvar=TRUE")
	if(!missing(by) && missing(data))
		stop("A data frame must be specified in order to use option by")
	if(missing(lead))
		stop("You must specify a lead")
	odata <- data
	ovar <- var	
	if(missing(by)){	
		var <- as.matrix(var)
		out <- as.matrix(var)
		for(i in (nrow(var)-lead):nrow(var)){
			out[i,] <- NA }
		for(i in 1:(nrow(var)-lead)){
			out[i,] <- var[i+lead,] }
		out <- as.vector(out) }
	if(!missing(by)){
		out.vector <- vector()
		for(z in 1:ncol(data)){
			#mess <- paste("data[,z] for z==",z," and lead ==", l, sep="")
			#print(mess)
			if(all(var == data[,z], na.rm=TRUE)){
				varnumber <- z}}
		pan <- unique(by)
		out <- vector()
		dat <- list()
		for (p in 1: length(pan)){
			dat[[p]] <- subset(data, by == pan[p])	}
		## This gets the two data sets needed
		for (p in 1: length(pan)){
			data <- dat[[p]]
			var <- as.matrix(data[,varnumber])
			out <- as.matrix(data[,varnumber])
			if(lead >= nrow(data)){
				for(i in 1:nrow(out)){
					out[i,] <- NA }}
			if(lead < nrow(data)){
				for(i in (nrow(var)-lead):nrow(var)){
					out[i,] <- NA }
				for(i in 1:(nrow(var)-lead)){
					out[i,] <- var[i+lead,] }}
			out.vector <- c(out.vector, out)}
	out <- out.vector}
	## The newvar code														
	if(newvar == TRUE){
		out <- as.data.frame(out)
		data <- cbind(odata, out)
		if(!missing(newname)){ 
			names(data)[[ncol(data)]] <- paste(newname) }
		if(missing(newname)){
			for(i in 1:(ncol(data) - 1)){
				if(all( (data[,i] == ovar) == TRUE, na.rm=TRUE )){
					varname <- names(data)[[i]] }}
			names(data)[[ncol(data)]] <- paste(varname, ".L", lead, sep="")	} 
		out <- data }	
	return(out) }


## Lag function
lag <- function(var, lag, by, data, newvar=FALSE, newname){
	if (missing(data))
		data <- environment(formula)
	if(newvar == TRUE && missing(data))
		stop("A data frame must be specified if newvar=TRUE")
	if(!missing(by) && missing(data))
		stop("A data frame must be specified in order to use option by")
	if(missing(lag))
		stop("You must specify a lag")
	if(missing(by)){
		tmp.data <- as.matrix(var)
		out <- as.matrix(var)
    lvar <- rep(NA, lag)
    lvar <- append(lvar, var)
    lag.matrix <- embed(var, lag + 1)
    lag.matrix <- lag.matrix[,-1]
		new.data <- as.matrix(lag.matrix) }
	if(!missing(by)){
    new.data<-data.frame()
    unique.by <-c(unique(by))
    data <- cbind(data, var)
    for (j in 1:length(unique(by))) {
      panel <- unique.by[j]
		  tmp.data <- data[which(by==panel),]
			lag.matrix <- matrix(NA,nrow=length(tmp.data[,1]),ncol=lag)
      lvar <- rep(NA,lag)
      lvar <- append(lvar, tmp.data$var)
      lag.matrix <- embed(lvar, lag + 1)
      lag.matrix <- lag.matrix[,-1]
      lag.matrix <- as.matrix(lag.matrix)
      new.data <- rbind(new.data, lag.matrix)}}   
	## The newvar code
	if(newvar == TRUE){
		if(!missing(newname)){
      colnames(new.data)<-paste("l",1:lag,".",newname,sep="")}
		if(missing(newname)){
      for(i in 1:(ncol(data) - 1)){
        if(all( (data[,i] == var) == TRUE, na.rm=TRUE )){
          varname <- names(data)[[i]] }} 
      colnames(new.data) <- paste("l",1:lag,".",varname, sep="")}	
    data<-cbind(data,new.data)}# end of newvar loop
  return(new.data) }


## Main lag variable. Some info in man page
alreg <- function(formula, data, model, method, lags, by, time.var, min.lag, max.lag, b0, B0, ml.type, mod.prob.method){

## SIMPLE TIME SERIES OR POOOLED ANALYSIS  
 if(method=="pooled"){
	# Ensure MCMC pack is installed
	check <- library()
	if(any(check$results[,"Package"] == "MCMCpack")) 
		require(MCMCpack)
	else
		stop("Please install MCMCpack using \n	install.packages(\"MCMCpack\")")
	# Direct data to the working environment if data object is missing from function call
	if (missing(data)){
		data <- environment(formula)
	}
	# Handle variable names 
	Terms <- terms(formula)                        
	vars <- as.list(attr(Terms, "variables"))
	vars[[2]] <- as.name("Intercept")
	varnames <- vars[-c(1)]
	final.vars <- vars[-c(1)]
	intercept.value <- attr(Terms, "intercept")
	# exception for no-intercept models
	if (intercept.value == 0){
		varnames <- varnames[-c(1)]
	}
	names(varnames) <- varnames
## Hardcoded prior? Do we really need this?
	wt <- rep(0.1, length(varnames)) 
	# Create and object for the rhs variables
	variables <- varnames
	variables <- variables[-c(1)]
	# set up triple matrices and create the shell of a results object
	for(q in 1:length(varnames)){
		if(min.lag == 0){
			triples <- matrix(nrow=(max.lag+1), ncol=3)}
		if(min.lag != 0){
			triples <- matrix(nrow=max.lag, ncol=3)}
		results <- list(triples, min.lag, max.lag, time)
		names(results) <- c("triples", "min.lag", "max.lag", "time")
		varnames[[q]] <- results
	}            
	print("...")
	print("This may take a few moments. Patience s.v.p.")
	
	
	
	
	
	mf <- model.frame(formula, data, na.action=NULL)
	# Creates a lag matrix for each of the variables we want lagged. New matrix is called "lagged[varname]". The lag mat is a N by (max.lag - min.lag) matrix
	# Then, bind this matrix to the mf object
	for(k in 1:length(lags)){                 
		name <- variables[[lags[[k]]]]
		lag.call <- paste("lagged",name," <- lag(data$",name,",",max.lag,",data=data,by=data$",by,")", sep="")
		eval(parse(text=lag.call))
		bind.call <- paste("mf <- as.data.frame(cbind(mf, lagged",name,"))", sep="")
		eval(parse(text=bind.call))
	}
	# NAs out all values for a row of the mf object with any missing values. This keeps the mf matrix conformable to the lagged[varname] mats and does not require alteration of the code to run the models. The NAs in the DV ensure that all models of all lags will have the same data after default casewise deletion
	for(p in 1:nrow(mf)){
		if(any(is.na(mf[p,]))){
			mf[p,] <- NA
		} 
	}



	# this loops through each of the lags computing each of the models
	for(l in min.lag:max.lag) {
		# This is added so each of the 4 variables are calculated simultaneously, ie each lagged variable has an equal lag
		for(d in 1:length(lags)){               
			name <- variables[[lags[[d]]]]
			# This line picks out the column (lag) from the matrix for each of the variables.
			# It replaces the unlagged variable in mf object with its appropriately lagged value from the lag object created in the k loop above.
			replace.call <- paste("mf$",name," <- lagged",name,"[,l]", sep="")
			eval(parse(text=replace.call))
			# The following lines implement the MCMCpack models for linear regression, logit, probit, and poisson respectively.
			if(model== "ls"){
				output.call <- paste("output",l," <- MCMCregress(formula, b0=b0, B0=B0, burnin = 1000, mcmc = 10000, marginal.likelihood= ml.type, data=mf)", sep = "")
				eval(parse(text=output.call))
			}
			if(model == "logit"){
				output.call <- paste("output",l," <- MCMClogit(formula,  b0=b0, B0=B0, burnin = 1000, mcmc = 10000, marginal.likelihood= \"Laplace\", data=mf)", sep = "")
				eval(parse(text=output.call))
			}
			if(model == "probit"){
				output.call <- paste("output",l," <- MCMCprobit(formula, b0=b0, B0=B0, burnin = 1000, mcmc = 10000, marginal.likelihood= ml.type, data=mf)", sep = "")
				eval(parse(text=output.call))
			}
			if(model == "poisson"){
				output.call <- paste("output",l," <- MCMCpoisson(formula,  b0=b0, B0=B0, burnin = 1000, mcmc = 10000, marginal.likelihood= \"Laplace\", data=mf)", sep = "")
				eval(parse(text=output.call))
			}		
			output.callback <- paste("output <- output",l, sep = "")
			eval(parse(text=output.callback))
			su <- output.summary(output)
		}
		# Fill in triple matrix for easy ropeladder plotting
		triples.mat <- cbind(su[,3], su[,1], su[,4])
		if(min.lag == 0){
			for(r in 1:(nrow(triples.mat))){
				varnames[[r]]$triples[l+1,] <- triples.mat[r,]
			}
		}
		if(min.lag != 0){
			if (model == "ls") {
				for(r in 1:(nrow(triples.mat)-1)){
					varnames[[r]]$triples[l,] <- triples.mat[r,]
				}
			}
			if (model == "probit") {
				for(r in 1:(nrow(triples.mat))){
					varnames[[r]]$triples[l,] <- triples.mat[r,]
				}
			}
		}
		notice <- paste("Finished Lag ",l,sep="")
		print(notice)
	}
	  
	# FUTUTE TO DO: Write more general code for BF. 
	
	if (mod.prob.method == "BF"){
		if (max.lag == 1) {
			BF<-BayesFactor(output1)
			postprob <- PostProbMod(BF)}
		if (max.lag == 2) {
			BF<-BayesFactor(output1, output2)
			postprob <- PostProbMod(BF)}
		if (max.lag == 3) {
			BF<-BayesFactor(output1, output2, output3)
			postprob <- PostProbMod(BF)}
		if (max.lag == 4) {
			BF<-BayesFactor(output1, output2, output3, output4)
			postprob <- PostProbMod(BF)}
		if (max.lag == 5) {
			BF<-BayesFactor(output1, output2, output3, output4, output5)
			postprob <- PostProbMod(BF)}
		if (max.lag == 6) {
			BF<-BayesFactor(output1, output2, output3, output4, output5, output6)
			postprob <- PostProbMod(BF)}
		if (max.lag == 7) {
			BF<-BayesFactor(output1, output2, output3, output4, output5, output6, output7)
			postprob <- PostProbMod(BF)}
		if (max.lag == 8) {
			BF<-BayesFactor(output1, output2, output3, output4, output5, output6, output7, output8)
			postprob <- PostProbMod(BF)}
		if (max.lag == 9) {
			BF<-BayesFactor(output1, output2, output3, output4, output5, output6, output7, output8, output9)
			postprob <- PostProbMod(BF)}
		if (max.lag == 10) {
			BF<-BayesFactor(output1, output2, output3, output4, output5, output6, output7, output8, output9, output10)
			postprob <- PostProbMod(BF)}
		if (max.lag == 11) {
			BF<-BayesFactor(output1, output2, output3, output4, output5, output6, output7, output8, output9, output10, output11)
			postprob <- PostProbMod(BF)}
		if (max.lag == 12) {
			BF<-BayesFactor(output1, output2, output3, output4, output5, output6, output7, output8, output9, output10, output11, output12)
			postprob <- PostProbMod(BF)}
		if (max.lag == 13) {
			BF<-BayesFactor(output1, output2, output3, output4, output5, output6, output7, output8, output9, output10, output11, output12, output13)
			postprob <- PostProbMod(BF)}
		if (max.lag == 14) {
			BF<-BayesFactor(output1, output2, output3, output4, output5, output6, output7, output8, output9, output10, output11, output12, output13, output14)
			postprob <- PostProbMod(BF)}
		if (max.lag == 15) {
			BF<-BayesFactor(output1, output2, output3, output4, output5, output6, output7, output8, output9, output10, output11, output12, output13, output14, output15)
			postprob <- PostProbMod(BF)}
	}

	if (mod.prob.method == "BIC"){
		BIC <- c()
		k <- length(final.vars)
		for (r in 1:max.lag){
			BIC.call<-paste("BIC[r]<- (-2)*(attr(output",r,", \"logmarglike\")) + k * log(length(attr(output",r,",\"y\")))", sep = "") 
			eval(parse(text=BIC.call))
			mBIC<- BIC - max(BIC)
			postprob <- exp(-0.5 * mBIC)/sum(exp(-0.5 * mBIC))
			postprob[is.na(postprob)] <- 1
			}
		}
			
	bma.mat<-matrix(nrow=length(variables)+1,ncol=max.lag)
	bma<-rep(NA, length(bma.mat[,1]))
	
	for(r in 1:(length(variables)+1)) {
		bma.mat[r,]<-varnames[[r]]$triples[,2]}
	
	bma.var.wt<-bma.mat
	
	for (m in 1:length(postprob)){
		bma.mat[,m]<-bma.mat[,m]*postprob[m]}
		
	for (r in 1:(length(variables)+1)){
		bma[r]<-sum(bma.mat[r,])}
	
	#pull se's/variance
	var.mat<-matrix(nrow=length(variables)+1,ncol=max.lag)
	
	for(r in 1:max.lag)	{
		var.call<-paste("vari<-(summary(output",r,")$statistics[,2])^2", sep = "")
		eval(parse(text=var.call))
		if (model == "ls") {
			var.mat[,r]<-vari[c(-length(vari))]}
		if (model == "probit") {
			var.mat[,r]<-vari}
		}
	
	#calculate mean posterior mean
	bma.mean<-c(rep(NA,length(bma.var.wt[,1])))
	for(r in 1:length(bma.var.wt[,r])){
		bma.mean[r]<-mean(bma.var.wt[r,])}
		
	variance<-c(rep(NA, length(bma)))
	for(r in 1:length(bma)){
		variance[r] <- sum(postprob*var.mat[r,])+ sum(postprob*((bma.var.wt[r,]-bma.mean[r])^2))
		}
	
	ub<-c()
	lb<-c()
	bma<-cbind(bma,sqrt(variance))
	if (model == "ls") {
		for (r in 1:length(bma[,1])){
			lb[r] <- bma[r,1]-1.96*bma[r,2]
			ub[r] <- bma[r,1]+1.96*bma[r,2]
			}}
			
	if (model == "probit") {
		alpha = 0.5
		for (r in 1:length(bma[,1])){
			lb[r] <- bma[r,1] - qnorm(1-alpha/2) * bma[r,2] 
			ub[r] <- bma[r,1] + qnorm(1-alpha/2) * bma[r,2] 
			}}
			
	bma<-cbind(bma, lb, ub)
	rownames(bma) <- final.vars
	colnames(bma) <- c("point estimate", "sd", "2.5%", "97.5%")
	
	varnames$min.lag <- min.lag
    varnames$max.lag <- max.lag
    varnames$model <- model
    varnames$formula <- formula
	varnames$bma <- bma
	varnames$postprob <- postprob
    class(varnames) <- "alpooled" 
	
	}
   
	
## PANEL ANALYSIS
   if(method=="bytime"){
    check <- library()
    if(any(check$results[,"Package"] == "MCMCpack")){
      require(MCMCpack)}
    else {
      stop("Please install MCMCpack using \n     install.packages(\"MCMCpack\")")}
    if (missing(data)){
      data <- environment(formula)}
    if(!is.vector(lags)){
      stop("lags must be a vector")}
    if(class(by) != "character"){
     stop("by must be of class character")}
    Terms <- terms(formula)
    vars <- as.list(attr(Terms, "variables"))
    vars[[2]] <- as.name("Intercept")
    varnames <- vars[-c(1)]
	final.vars <- vars[-c(1)]
    intercept.value <- attr(Terms, "intercept")
    if (intercept.value == 0){
      varnames <- varnames[-c(1)]}
    names(varnames) <- varnames
	variables <- varnames
    variables <- variables[-c(1)]
    time <- vector()
    time <- unique(time.var)
    time <- time[-c(1:max.lag)]
    tmp <- list()
                                
    pmean.mat <- sd.mat <- lb.mat <- ub.mat <- matrix(nrow=length(time), ncol=max.lag)
    results <- list(pmean.mat, sd.mat, lb.mat, ub.mat, min.lag, max.lag, time)
    names(results) <- c("pmean.mat", "sd.mat", "lb.mat", "ub.mat", "min.lag", "max.lag", "time")
    for(q in 1:length(varnames)){
      varnames[[q]] <- results}
    fo <- as.formula(print(formula))
    print("...")
    print("This may take a few moments. Patience s.v.p.")
    for(k in 1:length(lags)) {
      newdata <- data
      name <- variables[[lags[[k]]]]
      lag.notice <- paste("lagging ",name, sep="")
      print(lag.notice)
      lag.call <- paste("lagged",name," <- lag(newdata$",name,",",max.lag,",data=newdata,by=newdata$",by,")", sep="")
      eval(parse(text=lag.call))} # close lag loop
    for (l in min.lag:max.lag) { #for each of the lag lengths
      
	  for (d in 1:length(lags)) { #this sets it so that variables are lagged concurrently
        name <- variables[[lags[[d]]]]
        replace.call <- paste("newdata$",name," <- lagged",name,"[,l]", sep="")
        eval(parse(text=replace.call))}
        
	  for(t in 1:length(time)){   #for each of the subsetted ts
        tmp[[t]] <- subset(newdata, time.var == time[t])  
        time.notice <- paste("Time = ",t,sep="")
        print(time.notice)
		output.call <- paste("output",l," <- MCMCregress(formula, B0=B0, marginal.likelihood= ml.type, data=tmp[[t]])", sep = "")
		eval(parse(text=output.call))
        output.callback <- paste("output <- output",l, sep = "")
		eval(parse(text=output.callback))
		
		 
		if(model== "ls"){
			output.call <- paste("output",l," <- MCMCregress(formula, B0=wt, marginal.likelihood= ml.type, data=newdata)", sep = "")
			eval(parse(text=output.call))}
		if(model == "logit"){
			output.call <- paste("output",l," <- MCMClogit(formula, marginal.likelihood= \"Laplace\", data=newdata)", sep = "")
			eval(parse(text=output.call))}
		if(model == "probit"){
			output.call <- paste("output",l," <- MCMCprobit(formula, b0=b0, B0=B0, marginal.likelihood= ml.type, data=newdata)", sep = "")
			eval(parse(text=output.call))}
		if(model == "poisson"){
			output.call <- paste("output",l," <- MCMCpoisson(formula, B0=wt, marginal.likelihood= \"Laplace\", data=newdata)", sep = "")
			eval(parse(text=output.call))}
		
        output.callback <- paste("output <- output",l, sep = "")
		eval(parse(text=output.callback))
		  
		output<-output.summary(output)
		
        for(r in 1:nrow(uber)){
            varnames[[r]]$pmean.mat[t,l] <- su[r,1]
            varnames[[r]]$sd.mat[t,l] <- su[r,2]
            varnames[[r]]$lb.mat[t,l] <- su[r,3]
            varnames[[r]]$ub.mat[t,l] <- su[r,4]}	
	  }
      notice <- paste("Finished Lag ",l,sep="")   
      print(notice)	
	}
	
	
    varnames$time <- time
    varnames$min.lag <- min.lag
    varnames$max.lag <- max.lag
    varnames$model <- model
    varnames$formula <- formula
    class(varnames) <- "albytime" }
	
	return(varnames)
}

##PLOTTING FUNCTIONS

##ROPELADDER PLOTS
plot.alpooled <- function(object, titles, dims, margins, startvar, stopvar, ...){
	nvars <- length(object) - 4
	min.lag <- object$min.lag
	max.lag <- object$max.lag
	if(missing(titles)){
	titles <- names(object)
	tlength <- length(titles)
	titles <- titles[-c((tlength-3): tlength)]}
	par(mar=margins)
	par(mfrow=dims)
	for(i in startvar:stopvar){
		if(min.lag==0){
			plot(main=titles[i], las=1, xlab="", ylab="", (min.lag+1):(max.lag+1), object[[i]]$triples[,2], type="p",pch=19, col="black", cex=1.5, ylim=range(object[[i]]$triples), bty="n", xaxt="n")}
		if(min.lag!=0){
			plot(main=titles[i], las=1, xlab="", ylab="", min.lag:max.lag, object[[i]]$triples[,2], type="p",pch=19, col="black", cex=1.5, ylim=range(object[[i]]$triples), bty="n", xaxt="n")}
	for(j in min.lag:max.lag){
		if(min.lag==0){
			segments(j+1, object[[i]]$triples[j+1,1], j+1, object[[i]]$triples[j+1,3], lwd=2, col="gray50")
			points(j+1, object[[i]]$triples[j+1,2], cex=1.30, col="black", pch=19)}
		if(min.lag!=0){
			segments(j, object[[i]]$triples[j,1], j, object[[i]]$triples[j,3], lwd=2, col="gray50")
			points(j, object[[i]]$triples[j,2], cex=1.30, col="black", pch=19)}	}
	if(any(0 > object[[i]]$triples[,1]) && any(0 < object[[i]]$triples[,3])){abline(h=0, col="red")}   }}


## TOPOGRAPHICAL PLOTS
plot.albytime <- function(object, type=NULL, title, subtitle, coloring=NULL, colorby="t", scaleby=1, scale=TRUE, horizontal, vertical, ...){
	stackmat <- function(x){
		R <- nrow(x)
		C <- ncol(x)
		outmat <- x[-1, -1] + x[-1, -C] + x[-R, -1] + x[-R, -C]
		return(outmat) }
	C <- object$coef.mat
	SE <- object$se.mat
	T <- object$t.mat
	PRT <- object$prt.mat
	L <- object$lb.mat
	U <- object$ub.mat
	W <- object$wid.mat
	S <- object$sig.mat 
	max.lag <- object$max.lag
	time <- object$time
	if(type=="contour.coef"){
		x <- 1:nrow(C)
		y <- 1:ncol(C)
		filled.contour(x, y, C, color = terrain.colors,
			plot.title = title(main = print(title),
			xlab = "Time Period", ylab = "Lag"),
			plot.axes = { axis(1, seq(1, to=length(time), by = 1))
                  axis(2, seq(1, max.lag, by = 1)) },
			key.title = title(main="Size of \n Effect"),
			key.axes = axis(4, round(seq(min(C), max(C), by=((max(C)-min(C))/(10-1))), digits = 2) ))}
	if(type=="contour.t"){
		x <- 1:nrow(T)
		y <- 1:ncol(T)
		filled.contour(x, y, T, color = terrain.colors,
			plot.title = title(main = print(title),
			xlab = "Time Period", ylab = "Lag"),
			plot.axes = { axis(1, seq(1, length(time), by = 1))
                  axis(2, seq(1, max.lag, by = 1)) }, 
			key.title = title(main="t \n Value"), 
			key.axes = axis(4, round(seq(min(T), max(T), by=((max(T)-min(T))/(10-1))), digits = 2) ))}
	if(type=="contour.prt"){
		x <- 1:nrow(PRT)
		y <- 1:ncol(PRT)
		filled.contour(x, y, PRT, color = terrain.colors,
			plot.title = title(main = print(title),
			xlab = "Time Period", ylab = "Lag"),
			plot.axes = { axis(1, seq(1, length(time), by = 1))
                  axis(2, seq(1, max.lag, by = 1)) }, 
			key.title = title(main="p \n Value"), 
			key.axes = axis(4, round(seq(min(PRT), max(PRT), by=((max(PRT)-min(PRT))/(10-1))), digits = 2) ))}
	if(type=="significance"){ 
		z <- C/1000 
		t <- time 
		l <- (1:max.lag)
		fill <- matrix("green3", nr = nrow(z), nc = ncol(z))
		for(r in 1:nrow(S)){
			for(c in 1:ncol(S)){
				if(S[r,c] == 0){
					fill[r,c] <- "red3" }}}
		persp(t, l, z, theta = 0, phi = 90, scale = FALSE, col = fill, ltheta = -120, lphi = 15, ticktype = "detailed", xlab = "Time", ylab = "Lag", zlab = "", shade = 0.65, border = NA, box = TRUE)
		if(!missing(title)){mtext(print(title), side=3, line=2, cex=1.5)}
		if(!missing(subtitle)){mtext(print(subtitle), side=3, line=0.0, cex=1.5)}  }
	if(type=="persp"){
		z <- C * scaleby
		t <- min(time):max(time) 
		l <- (1:max.lag)
		if(colorby=="t"){CB <- T}
		if(colorby=="prt"){CB <- PRT}
		if(colorby=="se"){CB <- SE}
		if(colorby=="coef"){CB <- C}
		if(colorby=="ci.wid"){CB <- W}
		s.out <- stackmat(CB)
		if(is.null(coloring) || coloring=="heat"){fcol <- heat.colors(20)[cut(s.out, quantile(s.out, seq(0, 1, len = 21)), include.lowest = TRUE)]}
		if(!is.null(coloring) && coloring=="topo"){fcol <- topo.colors(20)[cut(s.out, quantile(s.out, seq(0, 1, len = 21)), include.lowest = TRUE)]}
		if(!is.null(coloring) && coloring=="cm"){fcol <- cm.colors(20)[cut(s.out, quantile(s.out, seq(0, 1, len = 21)), include.lowest = TRUE)]}
		if(!is.null(coloring) && coloring=="terrain"){fcol <- terrain.colors(20)[cut(s.out, quantile(s.out, seq(0, 1, len = 21)), include.lowest = TRUE)] }
		if(missing(horizontal)){horizontal = 100}
		if(missing(vertical)){vertical = 15}
		persp(t, l, z, scale = scale, col = fcol, theta=horizontal, phi=vertical, ticktype = "detailed", xlab = "", ylab = "Lag", zlab = "", shade = 0.65, border = NA, box = TRUE)
		if(!missing(title)){mtext(print(title), side=3, line=2, cex=1.5)}
		if(!missing(subtitle)){mtext(print(subtitle), side=3, line=0.0, cex=1.5)}		}
}


plot.almodel <- function(x, xlab = "", user.par = FALSE, ...) {
    if(class(x) == "alpooled" ){
    if (exists(paste("plot.alpooled"))){ UseMethod("plot.alpooled") } }
	if(class(x) == "alpanel"){
	if (exists(paste("plot.albytime"))){ UseMethod("plot.albytime") } }
    else
      stop(paste("class", class(x), "does not have an applicable plot method."))   }
	  
